Specific lieat of tlie Coulomb glass 

By A. Mobius 

Institute for Solid State and Materials Research Dresden, 
D-01171 Dresden, Germany, 

P. Thomas 

Department of Physics and Material Sciences Center, 
Philipps University, Renthof 5, D-35032 Marburg, Germany, 

J. Talamantes 
Department of Physics, California State University, 
Bakersfield, CA 93311, USA, 

and C.J. Adkins 
Cavendish Laboratory, 
Madingley Road, Cambridge CB3 OHE, UK 

April 30, 2001 



Abstract 

The specific heat c of the Coulomb glass is studied by numerical simu- 
lations. Both the lattice model with various strengths of disorder, and the 
random-position model are considered for the one- to three-dimensional 
cases. In order to extend the investigations down to very low temper- 
atures where the many-valley structure of the configuration space is of 
great importance we use a hybrid-Metropolis procedure. This algorithm 
bridges the gap between Metropolis simulation and analytical statistical 
mechanics. The analysis of the simulation results shows that the correla- 
tion length of the relevant processes is rather small, and that multi-particle 
processes yield an essential contribution to the specific heat in all cases 
except the one- dimensional random-position model. 



1 Introduction 

The Coulomb glass, a prominent model of disordered systems of localised par- 
ticles with long-range interaction, has been extensively studied for almost three 
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decades. It plays an important role as a semi-classical model for dilute im- 
purity bands when quantum interference as well as thermal excitation to the 
conduction band or from the valence band can be neglected. Various aspects of 
this model have been investigated since the pioneering works by PoUak (1970), 
Srinivasan (1970), Efros and Shklovskii (1975), and Kurosawa and Sugimoto 
(1975); for related reviews see Shklovskii and Efros (1984), Efros and Shklovskii 
(1985), Pollak and Ortuho (1985), and Ortuho et al. (2001). Nevertheless, this 
area of research has not lost its attraction, as is shown, for example, by recent 
studies of aging effects in electric conduction (Vaknin ct al. 2000) , of the possi- 
bility of a phase transition (Grannan and Yu 1993, Vojta and Schreiber 1994, 
Diaz-Sanchez et al. 2000b), and by combined investigations of dc conductivity 
and tunnel conductance (Sandow et al. 2001). 

An important and characteristic feature of the Coulomb glass is the many- 
valley structure of its configuration space, not only concerning single-particle 
hops, but also with respect to some classes of more complex excitations. In 
this respect it resembles the Ising spin glass and other disordered systems, see 
e.g. Mezard et al. (1987). This configuration-space structure is related to the 
existence of a huge number of so-called pseudo-ground states, as was observed 
in the first numerical Coulomb glass simulations (Baranovskii et al. 1976). The 
spectrum of those pseudo-ground states has recently been investigated in detail 
by Kogan (1998). Moreover, the related importance of many-particle excitations 
is directly illustrated by Fig. 3 in Schreiber and Tenelsen (1993). In particular, 
the huge number of local minima has been shown to cause a glass like relaxation 
behaviour (Schreiber et al. 1996, Perez-Garrido et al. 1999) and non-ergodic 
effects (Diaz-Sanchez et al. 2000a). 

However, under certain conditions, also "equilibrium" obscrvables can reflect 
the many-valley character of the configuration space. We focus here on the spe- 
cific heat c (heat capacity per site), in particular on its temperature dependence 
c(T) for various dimensions d. It has been studied numerically by several groups 
(Baranovskii ct al. 1982, Davics ct al. 1982, Davics ct al. 1984, Schreiber and 
Tenelsen 1993, Tenelsen and Schreiber 1994, Mobius and Pollak 1996, Mobius 
and Thomas 1997), mainly by means of the Metropolis procedure (Metropolis 
ct al. 1953). However, this approach is inappropriate in the low-temperature re- 
gion. One of the reasons for its failure is a direct consequence of the many-valley 
structure of the configuration space: the simulations usually get stuck in some 
local minimum. Consequently, other methods have to be used. Thus, almost 
complete sets of system states were obtained for the direct calculation of the low- 
T behaviour of c(T) by moans of different algorithms in Schreiber and Tenelsen 
(1993), Tenelsen and Schreiber (1994), Mobius and Pollak (1996), Mobius and 
Thomas (1997), compare also Mochena and Pollak (1991) and Perez-Garrido et 
al. (1997). Unfortunately, these specific heat studies have so far been restricted 
to the three-dimensional lattice model. 

To obtain an analytical description of c(T), Baranovskii et al. (1980) consid- 
ered the excitation from the ground state by independent single-particle hops. 
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In the case of randomly positioned sites, the adequacy of this approach was 
supported by numerical experiments (Baranovskii et al. 1982). However, for 
the three-dimensional lattice case, Mobius and PoUak (1996) showed (though 
only for one value of the disorder strength) that this description considerably 
underestimates c(T). 

The aim of our present work is to study the temperature dependence of the 
specific heat c(T) in the general case, considering in particular also that T region 
which is not accessible to usual Metropolis simulations. We have simulated one- 
to three-dimensional samples with both random and lattice structures, and, 
in the latter case, we have varied the strength of the disorder. To illuminate 
the puzzling situation described in the previous paragraph, we have sought 
for the conditions under which c{T) can be understood in terms of independent 
single-particle hops. For that purposes, the computer code which was developed 
in Mobius and Thomas (1997) for the three-dimensional lattice case has been 
rewritten and improved to handle arbitrary structures. 



2 Model 

All situations which are to be studied can be described by the Hamiltonian 

2^ \Xa-Xf3\ 
a a^P 

where the first sum represents an on-site disorder contribution and the second 
sum describes the Coulomb interaction, see e.g. Shklovskii and Efros (1984). 
The sums are over the N sites of the sample. These sites positioned at Xa 
are arranged either on a lattice (lattice model), or randomly (random-position 
model), the more realistic case. The occupation number e {1,0} indicates 
whether or not site a is occupied by one of the mobile charges - note that, in a 
donor band, Ua = I corresponds to a missing electron. In the lattice model, the 
Ea describe only the on-site disorder potential; they arc uniformly distributed 
between —B/2 and B/2 where the parameter B characterises the strength of 
the disorder. In this case, for elect roneutrality, background charges ~K are 
attached to each site where K = {ua). However, for an n-type impurity band 
emulated by means of the random-position model, the denote the potentials 
of the charged acceptors and we have K — 0. Elementary charge, spatial density 
of the Xa, and dielectric constant are all taken to be 1. 

In order to reduce finite-size effects and to exclude the influence of surfaces, 
wc have modified the interaction term in our simulations by imposing periodic 
boundary conditions and using the minimum image convention (Metropolis ct 
al. 1953, Davies et al. 1982). Thus our sample is an elementary cell in a (su- 
pcr)lattice, and the interaction energy of two sites is assumed to be determined 
by the shortest distance within this repeated representation. Moreover, we take 
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into account particle exchange with a reservoir. That means we consider the 
grand canonical potential h = H — nJ2a where the chemical potential n, if 
it is not already fixed by symmetry reasons, is determined by relaxing the closed 
system before starting the Metropolis part of the simulation. Considering the 
grand canonical ensemble rather than the canonical has several advantages. Cor- 
relation times arc reduced by emulating in this manner also long-range hops. 
Moreover, in particular at high T, the influence of finite-size effects is dimin- 
ished. However, one must remember that any T-dependence of ^ may influence 
c(T), as it does in metals. We have taken care to exclude such effects in the 
tests of our code. 

3 Algorithm 

In order to treat correlations completely, we have studied the specific heat of 
the Coulomb glass by numerical simulations. The standard procedure for this is 
the Metropolis algorithm (Metropolis ct al. 1953) which consists of substituting 
time averaging for ensemble averaging. For related surveys see e.g. Binder (1984) 
or Binder (1997), and references therein. However, in this approach, at low T 
and close to phase transitions, large difficulties arise from diverging correlation 
times. In the low-T case, these difficulties are reduced to a large extent by the 
hybrid algorithm developed in Mobius and Thomas (1997). This bridges the 
gap between Metropolis simulation and analytical statistical mechanics. In the 
present work, we use an improved version of this procedure. 

The difficulties mentioned above are tackled here in two ways following 
Mobius and Thomas (1997). First, in order to handle the problem of the great 
portion of unsuccessful attempts at low T, we do not perform random trials 
of state modifications at discrete time values as in the original Metropolis ap- 
proach. Rather, on a continuous time scale, wc focus on the state modifications 
which are really performed, though each such state investigation is clearly more 
elaborate than the trials in the standard procedure: To obtain the mean dwell 
time at the current state, the transition probabilities to all of its neighbours in 
the configuration space are calculated; the next state (in time) is picked ran- 
domly according to the ratios of the various transition probabilities, see Mobius 
and Thomas (1997). This approach resembles the n-fold way algorithm (Bortz 
et al. 1975). It becomes increasingly advantageous in comparison to the stan- 
dard Metropolis procedure as T is lowered. Moreover, problems which can arise 
from time discretization (Adam et al. 1999) are avoided in this way. 

Second, our hybrid algorithm takes advantage of the fact that thermody- 
namic averages are independent of the dynamics considered provided only that 
detailed balance is guaranteed. Thus, in order to reduce the problems arising 
from the huge number of local minima, it decomposes the configuration space 
into low- and high-energy parts, and utilises very fast auxiliary transitions ad- 
ditionally introduced between the states in the low-energy part. Therefore, the 



4 



sample simulated thermalises very rapidly with respect to the low-energy part 
of the configuration space as soon as its state enters this region: Due to the 
fast auxiliary transitions, the dwell time in this low-energy part as a whole is 
large in comparison to the dwell time in one of the states contained in it. Thus 
all quantities can be assumed to have the corresponding thermal averages of 
the low-energy subset of the configuration space while the sample stays in this 
configuration space region. Moreover, the escape probability from this subset 
does not depend on the entrance state. For details we have to refer to Mobius 
and Thomas (1997). For related approaches see Frantz et al. (1990) and Krauth 
and Pluchery (1994). 

In the present work, we have modified these ideas in three ways. We have 
used an improved algorithm for obtaining an almost complete set of low-energy 
system states. In Mobius and Thomas (1997), this subset was constructed by 
first sampling local minima, and then stepwise improving this subset by adding 
all the low-energy states found during the simulation, see also Schreiber and 
Tenelsen (1993). Here, we have extended the subset construction as introduced 
in Diaz-Sanchez et al. (2000a). That means, we first sample local minima; then 
we improve this backbone of the low-energy subset by thermal cycling (Mobius 
et al. 1997) (one could also compare this optimisation procedure with com- 
paction of a granular medium, such as ground coff'ee, by shaking). Finally, we 
search systematically the surroundings of all low-energy system states that have 
been found so far for further states (Mobius and Pollak 1996, Perez-Garrido et 
al. 1997). Also in our new code, we add the states visited during the Metropolis 
simulation to the low-energy states set if their energy is lower than the highest 
energy value in this subset. That this happens very seldom in practice, con- 
firms the efficiency of the three-stage low-energy subset construction procedure, 
described above. 

As has often been stressed by Mike Pollak, see e.g. the review Pollak and 
Ortufio (1985), cluster excitations are expected to play an essential role in the 
low-T dynamics of the Coulomb glass. In a direct way, one can include only 
small, compact clusters into the Metropolis dynamics. Unfortunately, clue to 
the long-range character of the interaction, the usual cluster algorithms as pro- 
posed in Swcndsen and Wang (1987) seem not to be applicable. Therefore, 
we use another approach here: Comparing repeatedly small rnimbers of states 
from the low-energy subset (mainly picked randomly), we search for groups of 
sites which change occupation only simultaneously, see Kawashima and Suzuki 
(1992). These clusters are tabulated and included into the simulation dynamics. 

Finally, in order to prevent the low-T improvements, described above, from 
slowing down too much the computer code in the high-T case, we have intro- 
duced a T dependent neighbourhood pruning: The number of the transitions 
included in the Metropolis dynamics is appropriately reduced with increasing 
temperature. 
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4 Numerical results 



A series of tests has been performed to check and optimise the new code. The 
quaUty of the sets of low-energy states was confirmed for various degrees of 
disorder by series of runs in which we treated the same sample repeatedly. 
Moreover, we simulated lattice samples without disorder and reproduced the 
ordered ground states. In these tests, when we considered half- filled three- 
dimensional samples of (2™)'* sites (m integer), we observed the NaCl structure 
to be the ground state only for odd m. The six-fold degenerate layered structure 
has slightly lower energies for arbitrary even m. At present, we cannot answer 
the question whether this strange feature of the long-range Coulomb interaction 
is only an artifact of the boundary conditions, or whether it has some physical 
meaning. However, it is remarkable in this context that layered structures 
have been observed in quench simulations of one-component plasmas (Ogata 
and Ichimaru 1989). Additionally, we have considered lattice samples without 
disorder where only a fraction 2"^ of the sites is occupied {p integer), studied 
analytically in Himt and PoUak (1986). For d = 3, p = 2, and m = 4, we 
obtained ground states with layered structure. 

We turn now to the discussion of the c(T) data obtained in the simula- 
tions. We have compared our results on the random-position model with those 
obtained previously by Baranovskii et ai. (1982): we observed very good agree- 
ment except for slightly larger peak values in our case. (In the present study, 
the smallest considered T value is lower than in Baranovskii et al. (1982) by 
a factor 16.) Furthermore, also the comparison with lattice data from Mobius 
and Thomas (1997) showed perfect agreement within the statistical accuracy. 

Finite-size effects can considerably complicate Coulomb glass simulations 
due to the long-range character of the interaction. Such problems are well-known 
from the study of the single-particle density of states where the sample size is 
related to a cut-off energy of the interaction spectrum and thus to a reliability 
limit in the energy dependence of this density of states (Baranovskii et al. 1979, 
Mobius et al. 1992). In consequence, the thermal smearing of the Coulomb 
gap is connected with a correlation length diverging as T ^ (Hmit 1990). 
However, those correlations and corresponding processes which are essential for 
the formation of the Coulomb gap need not also be the most important for the 
specific heat. Other state modifications rather than single-particle hops can 
be far more important for this quantity. This assertion is suggested by the 
numerical data from Mobius and Pollak (1996) showing that, for d = 3, c{T) 
data seem to be reliable down to T values far smaller than the cut-off energy of 
the interaction spectrum. 

To illuminate this point, we have plotted c{T) for various sample sizes in 
Fig. 1, considering the lattice model with fixed disorder strength with d = 1, 
2, and 3. For d = 1, the curves for N = 12, and 400 sites agree perfectly 
(within the statistical accuracy of our data) above as well as below the cut-off 
of the interaction energies, 2/A'', of the smaller sample. On the other hand, 
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considerable finite-size effects are present for iV = 3 over the range T < 1. Thus 
the specific heat is mainly determined by short-range single- or multi-particle 
hops here. There are large c(T) fluctuations in the smaller samples at the 
lowest T values (Tenelsen and Schreiber 1994, Mobius and PoUak 1996). These 
fluctuations seem to be "localised" in the big samples. 

In the case d = 2, significant finite-size effects arc present for TV = 9 and 
A'' = 25. However, the deviation between the data for samples of 100 and 900 
sites is very small, just at the accuracy limit. For d = 3, finite-size effects are 
clearly larger than in the previous cases. However, the results for samples of 
512 and 1000 sites agree well with each other. 

The convergence properties shown in Fig. 1 yield a rough estimate of the 
correlation radius: Zcorr = 3... 5. This length seems to increase "slowly" with 
decreasing T. 

Fig. 1 shows that the finite-size effects, which arise from the Coulomb in- 
teraction, rise with increasing d. To illuminate the interplay of interaction and 
dimension directly. Fig. 2 compares c(T) data for various d from lattice simu- 
lations with fixed disorder strength with the interaction-free case. In the latter 
situation, analytical calculations yield that c oc T for T — > and that c oc 
if T ^ oo, in agreement with the numerical data. Switching on the Coulomb 
interaction causes c to decrease at low T and to increase at high T, the more 
the larger d is. In the low-T region, c(T) decreases faster than linear as T ^ 0. 
Though c(T) cannot be approximated by a simple power law in the low-T re- 
gion studied, our data raises the question of whether c{T) might tend to a linear 
dependence again if T is sufficiently lowered. 

The influence of disorder strength is studied in Figs. 3 and 4, showing logj^g ^ 
and c, respectively, in dependence on logj^Q T . The observation that, for high T, 
increasing B reflects as some shift of c(T) towards higher T, which was reported 
for d = 3 in Mobius and Thomas (1997), is made here also for d = 2 and d=\. 
However, at low T , there is a weaker, but more complicated dependence on the 
disorder strength. In all cases, the random-position model yields broader c{T) 
curves than the lattice model. At low T, these data arc close to the results 
obtained for low disorder strength in the lattice model. 

A particular feature of c{T) is only seen clearly when c is plotted on a linear 
scale as in Fig. 4: In the lattice model, for B = 0.5 and d = 2 as well as d = 3, 
the c(T) curve is narrower and considerably higher than for the other B values. 
The question, whether or not this behaviour might be a precursor of a phase 
transition for still lower B (Grannan and Yu 1993, Vojta and Schreiber 1994, 
Dfaz-Sanchez et al. 2000b) cannot be answered yet. Corresponding simulations 
are extremely CPU-time consuming. 

Finally, we address the question of which processes determine c{T) as T — > 0. 
To clarify which contributions to the specific heat yield single-particle hops 
(SPH) and multi-particle hops, respectively, we have calculated the value csph, 
which one obtains when considering SPHs to be the only relevant, independent 
excitations of the ground state. The T-dependence of the quotient of csph and of 



7 



the "exact" c(T) data is represented in Fig. 5 for various d and B, for the lattice 
as well as for the random-position model. The SPH approximation obviously 
begins to fail when T is sufficiently increased for a significant number of sites 
to belong simultaneously to more than one of the excited SPHs. Therefore, 
Fig. 5 focuses on that region, where the corresponding overlap correction does 
not exceed 0.1c; quotient values for higher T, included in certain cases only, 
show that the breakdown of the above condition is correlated with a qualitative 
change of the quotient's T dependence. The maximum applicability T of the 
SPH approximation decreases with increasing sample size. Therefore, for d = I 
and 2, medium sample sizes are considered here. Nevertheless, quotient values 
could only be obtained for a T region which cannot be studied by standard 
Metropolis simulations, and, for the random-position model in the case d = 1, 
the applicability region of the SPH approximation includes only one data point. 

Fig. 5 shows that the SPH approximation, when applicable, clearly underes- 
timates c for all but one case; the random-position model for d = 1 is the only 
exception. In more detail, in the two- and three-dimensional lattice models, 
multi-particle contributions are particularly important for weak disorder. In 
the random-position model, the relevance of multi-particle excitations increases 
with the dimension. 

Thus, Fig. 5 implies the conclusion that, in agreement with Mike Pollak's 
view on the Coulomb glass model (PoUak and Ortufio 1985), multi-particle pro- 
cesses yield an essential contribution to the specific heat in the low-temperature 
limit. 
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Figure 1: Specific heat c (heat capacity per site) in dependence of temperature 
T for the lattice model with B = 1, K = 0.5. For the sake of clarity, c has 
been scaled by a dimension-dependent factor: qj = 1, 10~^, and 10~^ for d = 1, 

2, and 3, respectively. The symbols denote the the sample sizes as follows. +: 

3, 9, and 27, x: 12, 100, and 512, •: 400, 900, and 1000 for = 1, 2, and 3, 
respectively, and o: 25 and 125, for d = 2 and 3, respectively. All ensemble 
averages take into account such a large number of samples that the size of the 
la-error bounds does not exceed the symbol size. E.g., for the largest systems 
(•), we have considered 245, 92, and 50 samples for d = 1, 2, and 3, respectively. 



Figure 2: Dimensional dependence of c(T) for the lattice model with B = 2, 
K = 0.5. +: d=\,N = 400, x: o! = 2, TV = 900, and *: d = 3, = 1000. For 
comparison, results for interaction switched off, marked by •, are included. 



Figure 3: Influence of the nature and strength of disorder on c(T). For clarity, 
c has been scaled as in Fig. 1. +: B = 0.5, *: B = 2, and x: B = 8, all for the 
lattice model; o: random-position model. Error bars denoting the la deviation 
are only represented if their size exceeds that of the symbol. 



Figure 4: Influence of the nature and strength of disorder on c(T). +: B — 0.5, 
•: B = 1, *: B = 2, A: B = 4, and x: i? = 8, all for the lattice model; o: 
random-position model. The spline interpolations are guides to the eye only. 



Figure 5: Comparison of the specific heat value obtained in the single-particle 

hop approximation, cgpfj, with the c(T) values obtained when the correlations 
are taken into account completely. For the meaning of the symbols see Fig. 
4. The dashed curves represent csph data in cases, where the overlap of pairs 
cannot be neglected. For d = 1, 2 and 3, samples of 24, 400, and 1000 sites, 
respectively, were considered for the lattice model, and samples of 24, 400, and 
400 sites were simulated for the random-position model. 
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